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Abstract 

For odd square- free n > 1 the cyclotomic polynomial $„ (x) satisfies the identity of Gauss 

4<f„(x) = Al - (-l^-^nBl 

A similar identity of Aurifeuille, Le Lasseur and Lucas is 

<5> n {{-l) {n - 1)/2 x)=Cl-nxDl 

or, in the case that n is even and square-free, 

±$„/ 2 (-z 2 ) = Cl - nxDl, 

Here A n (x), . . . , D n (x) are polynomials with integer coefficients. We show how these coef- 
ficients can be computed by simple algorithms which require 0(n 2 ) arithmetic operations 
and work over the integers. We also give explicit formulae and generating functions for 
A n (x), . . . , D n (x), and illustrate the application to integer factorization with some numeri- 
cal examples. 
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1 Introduction 



For integer n > let $ n (x) denote the cyclotomic polynomial 



0<jr'<ra 
U,n)=l 

where £ is a primitive n-th root of unity. Clearly 



x n -l = l[* d (x), (2) 

d\n 

and the Mobius inversion formula [13] gives 

^n(x) = H(x d -ir^. (3) 

d\n 

Equation (1) is useful for theoretical purposes, but (3) is more convenient for computation 
as it leads to a simple algorithm for computing the coefficients of <& n (x) or evaluating Q n (x) at 
integer arguments using only integer arithmetic. If n is square-free the relations 

$ I ) - [ X ~ if n = 1; 

n[X) ~ I <5>n/p(x p )/<S>n/p(x), if p\n, p prime; [ ' 

give another convenient recursion for computing <& n (x). 

Although $ n (x) is irreducible over Z (see for example [29]), 3> n (x) may be reducible over 
certain quadratic fields. For example, 

4$ 5 (x) = (2x 2 + x + 2) 2 - 5x 2 , (5) 

so $5(x) has factors x 2 + ( 1=t= ^ ^ x + 1 whose coefficients are algebraic integers in Q[V5\. 
For odd square- free n > 1 the cyclotomic polynomial satisfies the identity 

4<M*) = ^-(-l) (n - 1)/ VB 2 . (6) 

Gauss [12] proved (6) for odd prime n; the generalization to other odd square-free n is due to 
Dirichlet [11]. Related identities of Aurifeuille and Le Lasseur [2] are 

• B ((-1)<«-M»x) = Cl - nxDl (7) 

for odd square- free n, and 

$> n ,2(-x 2 ) = C 2 n -nxD 2 n (8) 

for even square- free n > 2. For a proof, see Lucas [22] or Schinzel [26]. 

In (6-8), A n (x), . . . , D n {x) are polynomials with integer coefficients, and without loss of 
generality we can assume that A n (x)/2, B n (x), C n (x) and D n (x) are monic. In Section 3 we 
show how the coefficients of A n , . . . , D n can be computed by simple algorithms which require 
0(n 2 ) arithmetic operations and work entirely over the integers. 

In Section 1.1 we summarize our notation for future reference. Some numerical examples 
are given in Sections 1.2-1.3, and Newton's identities are discussed in Section 1.4. Then, in 
Section 2, we discuss the theoretical basis for the algorithms. The results for A n and B n are 
known (though perhaps forgotten) - they may be found in Dirichlet [11]. We present them in 
Section 2.2 for the sake of completeness and to aid the reader in understanding the results for 
C n and D n . 
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The algorithms are presented in Section 3. The algorithm (Algorithm D) for computing 
A n and B n is essentially due to Dirichlet [11], who illustrated it with some numerical examples 
but did not state it in general terms. The algorithm (Algorithm L) for computing C n and D n 
appears to be new. In Section 3.3 we comment briefly on Stevenhagen's algorithm [27] and 
compare it with Algorithm L. 

Finally, in Section 4 we give some explicit formulas for A n (x), . . . , D n {x). These may 
be regarded as generating functions if x is an indeterminate, or may be used to compute 
A n (x), . . . , D n (x) for given argument x. In the special case x = 1 the results for A n (l), B n (l) 
reduce to known formulas involving the class number of the quadratic field <3[\/±nJ. 

One application of cyclotomic polynomials is to the factorization of integers of the form 
a n ± b n : see for example [6, 7, 8, 9, 15, 16, 24, 26, 27]. If x = m 2 n for any integer m, then (7-8) 
are differences of squares, giving rational integer factors of x n ± 1. Examples may be found in 
Section 4.4. For the reader interested in integer factorization, our most significant results are 
Algorithm L of Section 3.2 and Theorem 3 of Section 4.4. 

1.1 Notation 

Unless qualified by "algebraic", the term "integer" means a rational integer, x usually denotes 
an indeterminate, occasionally a real or complex variable. 

jj,{n) denotes the Mobius function, 4>(n) denotes Euler's totient function, and (m, n) denotes 
the greatest common divisor of m and n. For definitions and properties of these functions, see 
for example [13]. Note that /i(l) = 4>(1) = 1. 

(m\n) denotes the Jacobi symbol 1 except that, as is usual for the Kronecker symbol 2 , (m\n) 
is defined as if (m,n) > 1. Thus, when specifying a condition such as (m\n) = 1 we may omit 
the condition (m,n) = 1. 

n denotes a positive integer (square-free from Section 2.2 on). For given n, we define integers 
n', s and s' as follows: 



valid for n = 2 (with C*2(x) = x + 1, D2{x) = 1). The Aurifeuillian factors of F n {x) are 






F+(x) = C n {x) + vW)„(z) 



and 




1 See, for example, Riesel [24]. To avoid ambiguity, we never write the Jacobi symbol as (^)- Note that m\n 
without parentheses means that m divides n. 
2 See, for example, Landau [18]. 
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From (10) we have F n {x) = F~(x)F+(x). We may write for one of F+, F~. 

We sometimes need to specify a particular complex square root. If m < then \fm means 
iy/\m\. 

d is usually the degree of a polynomial, while D is the discriminant of a quadratic form. For 
odd square-free n we always have D = sn, so D = 1 mod 4. 
Some additional notation is introduced in Section 1.4. 

1.2 Examples 

Taking n = 15, we have 

®is( x ) = !i = a: 8 - x 7 + x 5 - x 4 + x 3 - x + 1, 

(x & — l)(x d — 1) 

A 15 (x) = 2x 4 - x 3 - 4x 2 - x + 2, 

-Bi 5 (x) = x 3 - X, 

C 15 (x) = x 4 + 8x 3 + 13x 2 + 8x + 1, 

D 15 (x) = x 3 + 3x 2 + 3x + l, 

and the reader may easily verify that (6) and (7) are satisfied. As an example of (8), for n = 14 
we have 

F 14 (x) = $ 7 (-x 2 ) = = x 12 - x 10 + x 8 - x 6 + x 4 - x 2 + 1, 

x z + 1 

C 14 (x) = x 6 + 7x 5 + 3x 4 - 7x 3 + 3x 2 + 7x + 1, 

and 

D 14 (x) = x 5 + 2x 4 - x 3 - x 2 + 2x + 1. 

1.3 The identities of Beeger and Schinzel 

Taking n = 5 in (7) we obtain 

$ 5 (x) = (x 2 + 3x + l) 2 - 5x(x + l) 2 , 

so $5(x 2 ) has factors x 4 + 3x 2 + 1 ± \/5(x 3 + x) in Q[\/5]. Replacing x by x 3 we obtain factors 
of $ 5 (x 6 ). Now 

$i 5 (x 2 ) = 1> 5 (x 6 )/$ 5 (x 2 ), 
and by division (taking the factors with opposite signs of y/b) we obtain factors 

x 8 + 2x 6 + 3x 4 + 2x 2 + 1 ± \f5(x 7 + x 5 + x 3 + x) 

of $i5 (x 2 ). Thus 

$ 15 (x) = (x 4 + 2x 3 + 3x 2 + 2x + l) 2 - 5x(x 3 + x 2 + x + l) 2 . (11) 

This is not of the form (10) because it gives a factorization of $is(±x) over Q[\/ ±5] instead 
of Q[\/^15]. Instead, (11) is an example of the more general identities of Beeger [3] and 
Schinzel [26]. These identities can all be obtained in a similar manner from (10), so in the 
application to integer factorization they do not give any factors which could not be found from 
several applications of (10) and some greatest common divisor calculations. For this reason we 
have restricted our attention to identities of the form (6) and (10). 
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1.4 Newton's identities 

Let ^ 

P{ x ) = \{{ x -i j ) = Y,a 3 x d -i 

3=1 3=0 

be a polynomial of degree d with arbitrary roots £j and coefficients ao = 1, a±, . . . , a^- 
For k > 0, define 

Newton (1707) 3 showed how to express the elementary symmetric functions a±, a,2, ■ ■ ■ in terms 
of the sums of powers pi,p 2 , ■ ■ ■ 

We may find a\ , . . . , a& by solving a lower triangular linear system of special form [28] . 
Writing the solution explicitly in the form of a linear recurrence, we have 

k-i 

ka k = -^2 Pk~j a j (!3) 

3=0 

for k = l,...,d. An alternative expression for a k as a determinant may be obtained by applying 
Cramer's rule to the lower triangular system. However, for computational purposes (13) is more 
convenient. 

In Section 4 we use the following generating function [25] for (ao, ai, . . .): 

d I oo \ 

x d P(l/x) = ]T a 3 x j = exp - £ j (14) 
j=o V i =1 / 

Differentiating both sides of (14) and equating coefficients shows that (13) and (14) are formally 
equivalent. An independent proof of (14) is the following: for sufficiently small x we have 

d d oo oo / d \ oo 

in(^p(i/x)) = £ in(i - a kX ) = -£E ^*vj = - E E *J K'/j = - E l'r ,J i- 

k=l k=lj=l 3=1 \k=l ) j=l 

In all our applications of (14) the pj are bounded, so the infinite series converges for \x\ < 1. 

In the following, pj and aj are not fixed, but depend on the particular polynomial under 
consideration at the time. This should not cause any confusion. 

2 Theoretical Basis for the Algorithms 

Our idea is to compute sums of powers of certain roots of the polynomials occurring on the left 
side of (6) and (10), and then use Newton's identities in the form (13) to compute the coefficients 
of A n , . . . , D n . 

2.1 Cyclotomic polynomials 

First consider the computation of the coefficients of the cyclotomic polynomial & n (x) for n > 1. 
This is presented to illustrate a simple case of the technique; in practice it is more efficient to 
compute & n (x) from (3). 

3 See Turnbull [28], where the notation is used in place of our p^- It would be confusing to use Turnbull's 
notation because we have used s and s' for other purposes. Note that our pt are not generally prime numbers. 
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Let C be a primitive n-th root of unity. To apply Newton's identities we need to evaluate 



p k = E ( Jk ( 15 ) 



0<j<n 
C?»=l 

for k = 1,2, ... , 4>{n). This problem is well-known 4 . 
If n is prime the problem is easy: from 

1 — C n 

i+c+c 2 + ---+c n - 1 = TZ Y = ' 

we have p\ = —1. Moreover, for any k with (k,n) = 1, the map z i— > z k merely permutes 
{C, • • • ,C n_1 }, sop k =pi. 

Now consider the general case, n not necessarily prime. From (3) it is clear that p\ = /J,(n). 
Let gk = (k,n). If gk = 1 then the same argument as before shows that Pk = Pi- If <7fc > 1 then 
the sum (15) defining consists of <f>{n) / 4>{n / g^) copies of a sum of primitive (n/<7fc)-th roots 
of unity. Thus, the result is 

_ ri~/» W0 (16) 

Using (16) the coefficients ai, . . . ,o,M n ) °f c ^ri( a; ) may be evaluated from the recurrence (13). 

As an application of (14) and (16) we prove two Lemmas which give upper bounds on 
|<I > n (x)| and |F„(x)| for complex x outside the unit circle. Here F n (x) is the modified cyclotomic 
polynomial defined by (9). Lemma 2 is used in Section 4.4. 

Lemma 1 If \x\ > R > 1 then 

|* n (s)|<i2*Wexp(-^ 



Proof 

Let d = 4>(n). From (14) with x replaced by 1/x, we have 

$ n (x)/x d = exp ^- JpjaT J /jj , 

but from (16) we have 

\Pj\ < 9j < min(j,n), 

so 

\<f> n (x)/x d \ < exp ^S^'j = ex P {^~[ 
This completes the proof. □ 
Lemma 2 If n > 1 is square-free and \x\ > R > 1 i/ien 

|F n (x)| < J R^exp(-^l T ' 



4 If C = e 2 " 1 /™ then our pk is "Ramanujan's sum" c n (fc), in the notation of Ramanujan [23] or Chapter 26 of 
Davenport [10]. 
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Proof 

From the definition (9) of F n (x), 

J? _ J^( n )' if n is odd ; 

Q6g n \2<P(n/2), if n is even; 

so it is easy to see that 

deg F n = <t>{2n) 

in both cases. The bound on follows from Lemma 1 applied to $ n (±x) if n is odd, and 

from Lemma 1 applied to & n /2(— x 2 ) if n is even. □ 

2.2 The identity of Gauss 

From now on we assume that n > 1 is square-free. In this subsection we also assume that n is 
odd. Consider the polynomial 

G n (x)= I] ( X ~C J ) (17) 

0<j'<n 
(j\n)=l 

of degree 4>(n)/2, where £ = e 2 " K% l n is a primitive n-th root of unity. The particular choice of 
primitive root is only significant for the sign of the square root of sn appearing in the equations 
below. 

From Dirichlet [11], 

2G n (x) = A n (x)-^B n (x), (18) 



where A n and B n are as in (6), and s is defined in Section 1.1. Since n is odd, we have 
(-l|n) 
Define 



s = (-l| n ) = (_i)("-i)/2 



2G n (x) = A n (x) + yfiiiB n (x), (19) 
so Gauss's identity (6) may be written as 

$n(x) = G n (x)G n {x). (20) 

The sums of fc-th powers of roots of G n (x) are 



p k = E t Jk - (21) 

Let <7fc = (k,n). Then 



0<j<n 
(j\n)=l 



2pk = [ ^{n) + (k\n)^lE, if g k = 1; 
\n{n/g k )(l>(g k ), otherwise. 

The result (22) is essentially due to Dirichlet [11], but we sketch a proof. If g k = 1, then (22) 
follows from the discussion in Section 2.1 (where p k has a different meaning!) and the classical 
result that the Gaussian sum 

E owe' 

0<7<n 

is ^/sn. On the other hand, if g k > 1, observe that (gk, n /9k) = 1 because n is square-free. 
Thus, we can write the summation index j in (21) in the form j = jog k + ji(n/g k ), and (j|n) = 
(ji|5'fc)(jo|("/5 f fc))- Since jk = jogkk mod n, is independent of j±, and it follows that the 
sum (21) defining p k consists of (j}(g k )/2 copies of a complete sum of primitive (n/<7fc)-th roots 
of unity. Thus, (22) follows as in the proof of (16). 
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Although (22) has been written with two cases for the sake of clarity, our convention that 
(k\n) = if (k, n) > 1 implies that the expression 

2p k = fi(n/g k )(/)(g k ) + (k\n)y/sn (23) 

is valid in both cases. Similarly for G n (x), with the sign of \fsn in (23) reversed. 

Observe that p k G Q[\/^] is real if n = 1 mod 4, but p k £ QW~ n\ is complex if g k = 1 and 
n = 3 mod 4. 

Using (22), the coefficients of G n (x), and hence of A n (x) and B n (x), may be evaluated from 
the recurrence (13). Moreover, it is possible to perform the computation using only integer 
arithmetic. Details are given in Section 3. 



2.3 The identities of Aurifeuille, Le Lasseur and Lucas 

Here we assume that n > 1 is square-free, but not necessarily odd. Recall the definitions of n', 
s and s' from Section 1.1. 

Let C = e m / n ' be a primitive 2n'-th root of unity. The particular choice of primitive root is 
only significant for the sign of the square root in (29). Consider the polynomial 

L n (x)= U(x-e), (24) 

jeS„ 

where 

{j | < j < 2n>, (j, n') = 1, (j\n) = (-l)J}, if n = 1 mod 4; 

(25) 

{j I < j < 2n', (j, n') = 1, (n|j) = 1}, otherwise. 

Observe that L n (x) has degree 4>(n') = <p(2n). Also, j G 6* n iff 2n' —j G S'n, so the coefficients 
of L n (x) are real. In fact, from (29) below, they are in Qf-^/n]. We later use the fact that L n (x) 
is symmetric. 

Schinzel [26] essentially shows (with a different notation) that 

L n (x) = C n (x 2 ) - s'xV^D n (x 2 ) (26) 

where C n (x) and D n {x) are the polynomials of (10). Define 

L n (x) = L n (-x) = C n (x 2 ) + a'xy/nD n (x 2 ), (27) 

so after a change of variable (10) may be written as 

F n (x 2 ) = L n (x)L n (x). (28) 

Clearly F~(x) = L n (s' 'y/x) and F+(x) = L n (s'y/x). 

Let ^ = (k,n r ). The sums p& of fc-th powers of roots of L n (x) are 

_ J (n|/c)s'^/n, if k is odd; . . 

Pk ~ \n(n'/g' k )<f>{g' k )coB{(n- 1) fcvr/4) , if fc is even. 1 j 

Observe that the cosine in (29) is or ±1, and depends only on n mod 4 and k/2 mod 4. 
The proof of (29) is similar to that of (22), but tedious because of the number of cases to be 
considered. Thus, we omit the details. 

Using (29), the coefficients of L n (x), and hence of C n (x) and D n (x), may be evaluated from 
the recurrence (13). Details are given in Section 3. 
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3 Algorithms 

In this section we use the analytic results of Section 2 to derive efficient algorithms for computing 
the coefficients of the polynomials A n , . . . , D n . 

3.1 An algorithm for computing A n and B n 

Consider the computation of A n and B n for odd square-free n. Our notation is the same as in 
Section 2.2. Write 

d 

A n(x) = ^a jX d - j , 

3=0 
d 

3=0 

where d = <p(n)/2, qo = 2, (3q = 0, /3i = 1. 

Recall the definition (21) of p k . For k > we have, from (23), 



2pfc = qk + r k Vsn, 

where qt and r k are integers given by 

Qk = /J,{n/g k )(/)(g k ) (30) 

and 

r k = (k\n). (31) 
Using (13) and (18), we obtain the recurrences 

^ k-l 

a k = ^J2 ( snr k-jPj ~ Qk-jaj) (32) 

3=0 



and 



^ k-l 

Pk=y-Y^ ( r k-3 a 3 ~ Qk-jPj) (33) 
3=0 

for k = 1, 2, . . . , d. The algorithm is now clear: 
Algorithm D (for Dirichlet) 

1. Evaluate qk and r k for k = 1, . . . , d using (30-31). 

2. Set a <- 2 and /3 <- 0. 

3. Evaluate and for k = 1, . . . , d using (32-33). 

Comments on Algorithm D 

1. (32-33) should give exact integer results; in practice a sum not divisible by 2k is a symptom 
of integer overflow. 

2. The operation count can be reduced by a factor of close to four if advantage is taken of 
the following properties of A n and B n : 
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A n is anti-symmetric if its degree d = <f>(ri)/2 is odd, otherwise A n is symmetric (except 
for the trivial case A%(x) = 2x + 1). Thus, we may use 

oik = {-^) d a d -k 

H2k> d and n > 3. 

B n /x is antisymmetric if n is composite and n = 3 mod 4, otherwise B n /x is symmetric. 
Thus, we may use 

f fid-ki i n the symmetric case; 
fc I —fld-k, i n the anti-symmetric case. 

Using these properties, the recurrences (32-33) need only be applied for k < max(l, [d/2\). 
Example 

Consider the case n = 15 as in Section 1.2. We have s = —1, d = <^>(15)/2 = 4. Thus 
<7i = <72 = <?4 = ^(15)^(1) = 1, q 3 = m(5)0(3) = -2, 

n = (1|15) = 1, r 2 = (2|15) = (2|3)(2|5) = 1, r 3 = (3|15) = 0, r 4 = (4|15) = 1. 
33 1 34 1 r 3) and r4 are not required if we use symmetry. 

The initial conditions are «o = 2 and ft = 0. The recurrences (32-33) give 

a. x = (-15j"i/3 - 3i«o)/2 = -1, 

ft = (noo - giA))/2 = 1, 

a 2 = (-15r 2 ft - 15ri/3i - q 2 a - gi«i)/4 = -4, 

ft = (r 2 a + n«i - g 2 ft - gift)/4 = 0. 

Using symmetry of ^4 n (x) and anti-symmetry of B n (x)/x, or continuing with the recur- 
rences (32-33), we obtain 03 = a\ = —1, ft = —ft = —1, 04 = ao = 2, ft = —ft = 0. Thus 
^15 0*0 = 2x 4 — x 3 — Ax 2 — x + 2 and -Bis(x) = x 3 — x, as expected. 

3.2 An algorithm for computing C n and D n 

Consider the computation of C n and D n for square-free n > 1. Define n', s, s' and L n as in 
Section 2.3, and <i = <j){n')/2. Thus deg L n = 2d, deg C„ = d, and deg D n = d — 1. From (26) it 
is enough to compute the coefficients of L n (x). In order to work over the integers, we define 



Qk 



s'Pk/y/n, if k is odd; 
Pk, if & is even; 



where is the sum of fc-th powers of roots of L n (x). Thus, from (29), 

I (n\k), if A; is odd; 

9fe U(n7^)^)cos((n-l)A;7r/4), otherwise. ^ j 

If 

and 

z? ft (x)=x:^ 1 ^, 

J"=0 
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then, from (26), 



Ik = 0,2k 




7fc = ^E (nQ2k-2j-iSj ~ Q2k-2j7j) 



(35) 



and 



2k + 1 



1 



7fe + J! (<72fc+l-2j7i - Q.2k-2jSj) 
3=0 



k-1 



) 



(36) 



for k = 1,2,... 

We may use the fact that C n {x) and D n (x) are symmetric to reduce the number of times the 
recurrences (35-36) need to be applied. An algorithm which incorporates this refinement is: 

Algorithm L (for Lucas) 

1 . Evaluate q k for k = 1 , . . . , d using (34) . 

2. Set 70 <— 1 and <5 ^— 1. 

3. Evaluate j k for k = 1, . . . , [d/2\ and S k for A; = 1, . . . , |_(d - 1)/2J using (35-36). 

4. Evaluate 7fc for A; = [d/2\ + 1, . . . , d using jj, = jd-k- 

5. Evaluate 8^ for = |_(d + l)/2j , . . . , d — 1 using ^ = 



Consider the case n = 15 as in Section 1.2. We have n' = 2n = 30, s' = 1, d = </>(30)/2 = 4. 
Thus 



Using symmetry we obtain 73 = 71 = 8, 74 = 70 = 1, 62 = 5± = 3, and £3 = So = 1. Thus 
Ci 5 (x) = x 4 + 8x 3 + 13x 2 + 8x + 1 and £>i 5 (x) = x 3 + 3x 2 + 3x + 1, as expected. 



Example 



<Zi = (15|1) = 1, 
q 2 = A*(15)^(2) cos(7^) = -1, 

g 3 = (15|3) = 0, 
q 4 = M (15)^(2)cos(147r) = 1. 



The initial conditions are 70 = Sq = 1. The recurrences (35-36) give 



71 = (15gi<5 - ^270)72 = 8, 



Si = (71 + 9370 - <72<5o)/3 = 3, 
72 = (15<?3<5o + 15gi<5i - 9470 - <?27i)/ 4 = 13. 
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3.3 Stevenhagen's algorithm 



Stevenhagen [27] gives an algorithm for computing the polynomials C n (x) and D n {x). His 
algorithm depends on the application of the Euclidean algorithm to two polynomials with integer 
coefficients and degree 0(n). C n {x) and D n (x) may be computed as soon as a polynomial of 
degree < <p(n)/2 is generated by the Euclidean algorithm. Thus, the algorithm requires 0(n 2 ) 
arithmetic operations, the same order 5 as our Algorithm L. 

Unfortunately, Stevenhagen's algorithm suffers from a well-known problem of the Euclidean 
algorithm [14] - although the initial and final polynomials have small integer coefficients, the 
intermediate results grow exponentially large. When implemented in 32-bit integer arithmetic 
we found that Stevenhagen's algorithm failed due to integer overflow for n = 35. 

Algorithm L does not suffer from this problem. It is easy to see from the recurrences (35- 
36) that intermediate results can grow only slightly larger than the final coefficients 7& and 5k- 
A straightforward implementation of Algorithm L can compute C n and D n for all square-free 
n < 180 without encountering integer overflow in 32-bit arithmetic. When it does eventually 
occur, overflow is easily detected because the division by 2k in (35) or by 2k + 1 in (36) gives a 
non- integer result. 



4 Explicit expressions for A n , . . . , D n 

We now use (14) to give generating functions for the coefficients of A n , . . . , D n . The generating 
functions can be used to evaluate the coefficients of A n (x), . . . ,D n (x) in O(nlogn) arithmetic 
operations, via the fast power series algorithms of Section 5 of Brent and Kung [5]. Also, where 
the generating functions converge, they give explicit formulas which can be used to compute 
A n (x), . . . , D n (x) at particular arguments x. However, it is often more efficient to compute the 
coefficients of the polynomials by the algorithms of Section 3 and then evaluate the polynomials 
by Horner's rule. 

The generating functions may be written in terms of certain analytic functions f n and g n , 
which we now define. 

4.1 The analytic functions f n and g n 

For odd square- free n > 1 and |x| < 1, define 



/„(*) = £(j|n)y. (37) 

j=l J 

Similarly, for square- free n > 1 and |x| < 1, define 

g n { x ) = Y J {n\2j + l)^—_ -. (38) 

Observe that g n {x) is an odd function, so g n (— x) = —g n (x)- 

It follows from (46) and (57) below that exp(y / in/ n (x)) and exp(2y / n<7 n (a;)) are rational 
functions with zeros and poles at certain roots of unity. From these representations it follows 
that the analytic continuations outside the unit circle are given by 

fn{x) = f n (l/x), ifn = lmod4; . . 

/„(x) + /„(l/x) = 2/ n (l), ifn = 3mod4; { ™> 



5 The complexity of both algorithms can be reduced to 0(n(log n) 2 ) arithmetic operations by standard "divide 
and conquer" techniques [1, 4], but this is not of practical significance. 
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and 

g n (x) = g n {\/x). (40) 

The functions f n {x) and g n (x) are closely related. For example, taking the odd terms in the 
sum (37) and using the law of quadratic reciprocity, we obtain 

fnix) - fn(-x) = 2g n (xyfs)/yfs. 

Such identities are a consequence of relationships between the polynomials G n (x) and L n (x). 

/ ra (l) is related to the class number h(D) of the quadratic field QfV^D] with discriminant 
D = sn. In the notation of Davenport [10], / n (l) = L-i(l) = £(1) = L(l,x), where x(j) = U\ n ) 
is the real, nonprincipal Dirichlet character appearing in (37). Known results [10, 18, 19, 30] in 
the case n = 3 mod 4 (so D = —n) are 

"7T v-^ . ., . . 7T .c-^ , ., . Z7T 



/..(i) = ^ EWI-M = (2 _ (2 | n)) ^ g 01") - ^"(-») > o. . n 



Here 

n = 3; 

n = 3 mod 4, n > 3 



J 6, if 



is the number of roots of unity in Q[y/D]. Since h(—n) is an integer and /i(— 3) = 1, we have 

exp(2V^n(l)) = {5" 1 + ^ )/2 ' ^ n = 3 > AA (42) 
(1, if n = 3 mod 4, n > 3. 

In the case n = 1 mod 4, we have D = n, and 

/ n (l) = ^fc(n). (43) 



Here e is the "fundamental unit", i.e. e = (\u\ + y/n\v\)/2, where (u,v) is a minimal nontrivial 
solution of u 2 — nv 2 = 4. For example, if n = 5 then e = (3 + y/E)/2, h(5) = 1 and /s(l) = 
(In e)/ v / 5 = 0.4304... 

4.2 The polynomials A n and £> n 

Let n > 3 be odd and square-free. We exclude n = 3 to avoid the special case in (42), but the 
results apply with minor modifications when n = 3. Let s, G n , G n be as in Section 2.2, and 
d = 4>(n)/2. Recall that 

Gn(x)= I] (^-C J ) (44) 

0<j<n 
(j\n)=l 



and 



From (14) and (23), we have 



Gn(x)= [] ( X ~C J )- ( 45 ) 



0<j<n 
(j\n)=-l 



G n (l/x)/G n (l/x) = exp( v / ^/n(x)). (46) 

Also, from (44), 

(-x) d G n (l/x) = J] (C 3 x-1)= I] C j (*-C j ), (47) 

0<j<n 0<j<n 
(j\n)=l (j\n)=l 
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so 

G n (l/x)/G n (l/x) = C l[(x-C j )^ n \ (48) 

where 

n-1 

^=E0»J- (49) 

3=1 

If n = 1 mod 4 then by grouping the terms for j and n — j (j < n/2) in (49) we have n\a. If 
n = 3 mod 4 then, from (41), we have a = —nh(—n) so again n\a. Thus, in both cases ( a = 1, 
and from (48) we have 

r (1M/G (Mr) - f G n(x)/G n (x), ifn=lmod4; ( , 

Gn{1/X)/Gn{1/X) -\G n (x)/G n (x), ifn = 3mod4. (5 ° j 

It follows from (46) that 

G n (x)/G n (x) = exp(sy/snf n (x)). (51) 

We see from (46) or (51) that, as claimed above, exp( v / in/ n (x)) is a rational function. It has 
zeros at = —s; and poles at ( 3 , (j\n) = +s. From (20) and (51), taking a square root, 

we obtain 



G n (x) = J^{x) exp -^fHlf n (x) (52) 



If (52) is interpreted as a generating function for G n (x) then >/<J> n (x) and f n {x) should be 
interpreted as power series in x, and the correct sign of the square root is positive. On the 
other hand, if (52) is regarded as an exact expression for G n (x), then the sign of the square root 
is positive for real x, because G n (x) and <& n {x) have no real roots, and the exponential never 
vanishes, so a change in sign would contradict the continuity of G n (x). An extension of this 
argument shows that the same branch of the square root must be taken in any simply-connected, 
closed region which does not contain any of the zeros of <3? n (x). (We omit similar comments 
below.) 

From (52) we easily deduce the corresponding expressions for A n (x) = G n (x) + G n (x) and 
B n (x) = {G n {x) — G ra (x))/y / sn. We state the results as a Theorem: 

Theorem 1 For odd, square- free n > 3, the polynomials A n (x) and B n (x) occurring in Gauss's 
identity (6) are 

A n (x) = 2 v /<D n (x)cosh ^/nWj (53) 

and 



B n (x) = 2\ *^ sinh ( ^f n (x)) . (54) 



sn 



Remark 

If n = 3 mod 4, so s = — 1, then it is natural to replace cosh(iz) by cos(z) in (53) and sinh(iz) 
by isin(z) in (54), giving 



A n (x) = 2 x /$„(x) cos ( ^fn(x) (55) 



and 



B n (x) = 2\ ^1 sin f & n {x) ) . (56) 



n 
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Example 

Consider the case n = 15. We expand the right side of (55) as a power series in x, keeping 
enough terms to find A±^(x) without using symmetry. From Section 1.2 we have 

$ 15 (x) = 1 - x + x 3 - x 4 + x 5 - x 7 + x 8 , 



so 



x x 2 7x 3 37x 4 



rr ~~ x x (x" 
\/$15(x) = 1 1 



Also, 



2 8 16 128 



2 4 7 

X X X 



+ 



SO 



and 



/i 5 (x) = x + y + — - — + ••• 
/ , , A , 15x 2 15x 3 15x 4 

cos — /«(*) =i--g r + + 



2 v /$i 5 (x)cos ^^-/i 5 (x)j = 2 - x - 4x 2 - x 3 + 2x 4 + • • • , 

which is Ai§{x) + 0(x 5 ) as expected. We can ignore the 0(x 5 ) term (which in fact vanishes) 
since we know that deg (x) = </>(15)/2 = 4. 

A reader who attempts similar computations for larger n will soon be convinced that Algo- 
rithm D of Section 3.1 is more convenient, if only because all intermediate results are integers 
and there is an easy check on the accuracy of the inner product accumulations. 

Using (41-43) and the fact (an immediate consequence of (4)) that 



<Mi) 



n, if n is prime; 

1, for other square- free n > 1; 



we can verify that (46) and (51-56) give sensible results in the case x = 1. For example, we 
have 15) = 2, so /i 5 (l) = 27r/v^l5 and cos(^/i 5 (l)) = cos(tt) = -1. Thus (55) gives the 
correct value —2 of ^15(1). 

4.3 The polynomials C n and D n 

We now consider the analogues of (53-54) for the Lucas polynomials C n and D n . The argument 
is similar to that of Section 4.2, but simpler because the polynomial L n (x) is symmetric, which 
leads to the simple functional equation (40) for g n (x). 

Assume that n > 1 is square-free, and adopt the notation of Section 2.3. Using (29) in the 
generating function (14), we have 

L n (l/x)/L n (l/x) = exp (2s'^9n(x)) . (57) 

This shows that exp(2- v /ng r „(x)) is a rational function. From 

L n (l/x)/L n (l/x) = L n (x)/L n (x) 

we deduce the functional equation (40), which gives the analytic continuation of g n (x) outside 
the unit circle. We may also write (57) more simply as 

L n (x)/L n (x) = exp (2s' y/ng n (x)) . (58) 

From (28) and (58), taking a square root, we obtain 

L„(x) = ^F n (x 2 )exp (-s'^fng n (x)) . (59) 
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Theorem 2 Let n > 1 be square-free. The Aurifeuillian factors F^{x) = C n {x) ± y / nxZ? n (x) 
of F n (x) are given by 

F±{x) = y / F n (x) exp (±Vn-g n (Vx~)) . (60) 

Also, 

C n{x) = \J F n (x) cosh (y/ng n (y/x)) (61) 

and 

Dn(x) = J ^ sinh (Vn-g n (Vx~)) . (62) 
V nx 



Proof 

Recall that the Aurifeuillian factors F^{x) = C n {x) ± y/nxD n {x) of F n (x) are L^iy^). 
Thus (60) follows from (59). Since 

L n (x) + L n (-x) = 2C n (x 2 ) 

and 

L n (x) - L n (-x) = -2s'x^/nD n (x 2 ), 
we easily deduce (61-62). □ 

4.4 Application to integer factorization 

In this section we illustrate how the results of Sections 3.2 and 4.3 can be used to obtain factors 
of integers of the form a n ±b n . Our examples are for illustrative purposes, so are small enough 
to be verified by hand. Many larger examples can be found in [6, 7]. 

As usual, n > 1 is a square- free integer. Recall the definition of F n {x) from Section 1.1. 
Note that the polynomial F n (x) is a factor of x n ± 1 (where the sign is "— " if n = 1 mod 4, and 
"+" otherwise). 

If x has the form m 2 n, where m is a positive integer, then ^fnx = mn is an integer, and 
the Aurifeuillian factors F^{x) = C n (x) ± mnD n (x) give integer factors of F n (x), and hence of 
x n ± 1 = m? n n n ± 1. For example, if m = n k , we obtain factors of i7,( 2fc + 1 ) n ± 1. 

More generally, if m = p/q is rational, we obtain rational factors of x n ± 1 = p 2n q~ 2n n n ± 1, 
and thus integer factors of p 2n n n ± q 2n . We consider one example later, but for the moment we 
continue to assume that m is an integer. 

Before giving numerical examples, we state explicitly how the results of Section 4.3 can be 
used to compute 



m 2 n 



Theorem 3 Let m, n be positive integers, n > 1 square-free, x = m 2 n, and A = 0(2n)/2 . 
Then the Aurifeuillian factors F^{x) of F n {x) are given by 

F~(x) = [F+ 1/2J 

and 

F+(x) = F n (x)/F-(x), 

where 

1 ^ {n | 2j + 1) 



F = ^F n (x)exp\--J2 



mf^ Q {2j + l)xi 
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Proof 

From (60), using the functional equation (40) and the power series (38) for g n (l/^/x), we 
have 

F -(x) = Jf~(x) exp 1--Y ( ," I - 2J '| 1) 

Thus, we only have to show that the error incurred by truncating the power series after A terms 
is less than 1/2 in absolute value, i.e. that \F — F~\ < 1/2. 

If x < 5 we must have m = 1, n = 2 or 3, x = n, A = 1. In both cases F~{n) = 1 and it is 
easy to verify that 1/2 < F < 3/2. Thus, from now on we assume that x > 5. 

Let 

l_^(n|2j + l) 



Since m > 1, A > 1, x > 5, and the Jacobi symbol is or ±1, we have 

-A 



X~ 3 (\ 1 1 

l^g(2jTl)^(3 + 5T5 + 7^ + 



x 



so 

\t\ < 5 ^y/E In ^ 1 ~ *J X ~ X < 0.3802a;- A . (63) 

In particular, \t\ < 0.3802/x < 0.08, so 

/ x , 1*1 0.3802x- A \ 

exp (t) - 1 < ^— < < 0Ax- x . (64) 

1 FW 1 X — |*|/2 0.96 v ; 

Now F = F~ exp(t), so 

\F-F-\ = F-\exp{t)-\\. (65) 
Applying Lemma 2 of Section 2.1 with R = 5 gives 

F n (x) < exp(l/4)x^ 2n ), 

but (f)(2n) = 2A and F~ is the smaller factor of F n , so 

(x) < y/F n (x) < exp(l/8)x A . (66) 
From (64-66) we finally obtain the bound 

\F-F~\ < 0.4exp(l/8) < 0.5, 

which completes the proof. □ 

Since the argument of the exponential in Theorem 3 is —\/m+0(\/ri) as n — > oo, we have the 
following result, which sheds some light on the ratio of the Aurifeuillian factors. (The Corollary 
strictly follows from Theorem 3 only if m is an integer, but from the proof of Theorem 3 it is 
clear that the Corollary is also valid for rational m.) 

Corollary 1 Let x = m 2 n where m > is rational n > 1 is an integer. Consider m fixed as 
n — > oo through square-free values. Then the Aurifeuillian factors F^{x) of F n (x) satisfy 

F+{x)/F-{x) = exp(2/m) + 0(l/n). 
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Examples 

1. To start with a simple example, consider n = 2, m = 2, so x = m?n = 8 and A = 1. In 
Theorem 3 we have F 2 (x) = x 2 + 1 = 65, 

F = y / F 2 (x)exp(-l/m) = >/65exp(-l/2) = 4.89 ■■ ■ 

and rounding to the nearest integer gives the factor 5 of 65. 

2. A similar but less trivial example is n = 2, m = 2 5 , x = m 2 n = 2 11 . Here F 2 (x) = 
x 2 + 1 = 2 22 + 1 and 

F = yf F 2 (x) exp(-l/m) = V2 22 + lexp(-2~ 5 ) = 1984.98- ■■ 

which on rounding gives the factor F2 = 1985 of 2 22 + 1. By division we find F^~ = 2113, so 
the complete factorization is 2 22 + 1 = 5 • 397 -2113 . 

3. Now consider n = 5, m = 3, so x = m 2 n = 45 and A = </>(10)/2 = 2. In this case 

F 5 (x) = $ 5 (x) = (x 5 - l)/(x - 1) = 4193821, 

F = JF 5 (x)exp (-— + ^3-) = V4193821 exp(- 134/405) = 1470.99924 • • • 

and rounding to the nearest integer gives the factor 1471 of F$(x). By division we obtain the 
other factor 2851. In this example, but not in general, the Aurifeuillian factors are prime. 

4. Now consider a composite n, say n = 15. To keep the arithmetic easy we take m = 1, so 
x = 15 and 

F 15 (x) = $i 5 (-x) = x 8 + x 7 - x 5 - x 4 - x 3 + x + 1 = 2732936641. 
The example in Section 3.2 shows how we can use Algorithm L to compute 

C 15 (x) = x 4 + 8x 3 + 13x 2 + 8x + 1 

and 

D 15 (x) = x 3 + 3x 2 + 3x + l. 
Evaluating the polynomials, we find 

Ci 5 (15) = 80671, £>i 5 (15) = 4096, 

so the Aurifeuillian factors of Fi 5 (15) are 80671 ± 15 • 4096. This gives 19231 and 142111, 
which can easily be verified to be factors of 15 15 + 1. The "algebraic" factors (15 3 + 1)/16 and 
(15 5 + 1)/16 allow us to complete the factorization: 

15 15 + 1 = 2 4 • 31 • 211 • 1531 • 19231 • 142111 

Alternatively, instead of evaluating Cis(15) and -Dis(15), we can use Theorem 3. We have 
A = 0(15)/2 = 4. Since (15|3) = (15|5) = and (15|7) = 1, the computation gives 

F = V /Fi 5 (x)exp ^-1 - ^3^) = 19231.00217- ■■ 

and rounding to the nearest integer gives the factor 19231 of Fi5(15). 



18 



5. To conclude, we give an example where m = p/q is rational but not an integer. Consider 
n = 7, p = 2, q = 5, so x = p 2 n/q 2 = 28/25 and y/nx = pn/q = 14/5. We have 

j-i / \ ^ , x x 7 + 1 28 7 + 25 7 
F 7 (x) = *,(-*) = — = -^-^ 

Theorem 3 is not applicable. Because x is close to 1, the series for g 7 (\jyfx) converges rather 
slowly, and we need to take at least 35 terms to obtain sufficient accuracy. However, using 
Algorithm L we easily find that 

C 7 (x) = x 3 + 3x 2 + 3x + 1, D 7 (x) = x 2 + x + 1, 

so Horner's rule gives 

C 7 (x) = 148877/5 6 , D 7 (x) = 2109/5 4 , 

and C 7 {x) ± y/nxD 7 {x) gives 296507/5 6 and 1247/5 6 . Thus, we have obtained factors 296507 
and 1247 of (25 7 + 28 7 )/53. The larger factor is prime, so it is easy to deduce the complete 
factorization 

25 7 + 28 7 = 29 • 43 • 53 • 296507 
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